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Introduction 


The Monte-Carlo program BSIMUL was designed to simulate the interactions between 
gamma-rays and detectors. BSIMUL follows the path of a gamma-ray through a detector 
and can also be used to simulate spectra generated by a detector. Several improvements 
have been made to the program within the past year. The menu has been modified to allow 
the operator to use an input data file containing the simulation parameters. This enables the 
operator to set up a series of runs easily. Major portions of the program were rewritten to 
expand the possible detector configurations, and include tapered cylinders. Isotope and 
absorption coefficient libraries have been created, and the coefficients of compounds not in 
the coefficient library can be calculated if the density is known. 

The photopeaks generated by BSIMUL are one channel wide, sufficient for peak-to- 
Compton ratio calculations, but not sufficient for predicting the true shape of a spectrum. 
Standard peak shapes are available for several detector materials and can be used to add a 
true peak shape to BSIMUL's output. BSIMUL now writes the files necessary to perform 
an efficiency calibration with the spectrum fitting program ROBFIT. 

This report details the progress made in simulating gamma-ray spectra for a variety of user- 
specified detector designs. Several simulations were performed with the NEAR detector 
configuration. This report also includes a short discussion of work done in the related areas 
of spectral analysis and shield design in the performance of contract F08650-94-C-0006, for 
the Advanced Research Projects Agency and the U.S. Air Force. The work on shield 
design is presented here as an illustration of the capabilities of the program. 
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The Simulation Program 


A random number generator is used to simulate the behavior of gamma-rays as they pass 
through a detector. Absorption coefficients of the materials making up the detector, shield 
and collimator are used to calculate the reaction probability for: 1) photoabsorption; 2) 
Compton absorption; 3) pair production; and 4) Bremsstrahlung. The direction that the 
gamma-ray takes after a Compton interaction is calculated using the random number 
generator in conjunction with 


cos 0=1 



il 



[Ref 1], while the energy of the gamma-ray is calculated using the Klein- Nishina formula 
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BSIMUL was originally written with the assumption that the detector and shield would be 
cylindrical ly symmetric, and have constant radii. Removing the latter assumption created a 
substantial modification of the program, i.e., the subroutine used to calculate the time to a 
boundary between materials has been rewritten to enable the operator to model detectors 
with tapered parts. Rectangular parts can also be modeled by a slight modification to the 
collimator subroutines. 

Since the program was originally written to model germanium detectors, it was assumed 
that there would be a cold finger, dead layer, aluminum case, and possibly a shield. The 
detector material was assumed to be germanium, and changing it required editing the 
program. One array was used to specify both material and part, so that placing the shield 
inside the case, for instance, also required the program to be edited. The new version of 
the code separates the material from the part, and gives the operator complete freedom to 
use any design that contains fewer than 20 parts. One possible detector configuration is 
shown in Figure 1 . The collimator subroutines that were discussed in the last annual report 
have now been fully tested. 

In preparing to use the program with a graphical interface, the menu system has been 
replaced by an input file, shown in Table I. The first several lines (Table I-A) contain the 
number of events to be simulated, the isotopes and their relative strengths, and the source 
location. If the operator wants to simulate the spectrum of a source with a known strength, 
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it may be entered in nCi, jaCi, 
or mCi. The program will 
calculate the number of 
gamma-rays that will be 
emitted in the direction of the 
detector. A typical spectrum 
contains between one million 
and five million gamma-rays, 
only a few of which deposit 
energy in the detector. A 
spectrum of one hundred 
thousand gamma-rays is 
useful only for studying the 
photopeaks. The energy 
entered in line three is used 
only if a single energy 
spectrum is desired. If 
isotopes are specified, the 
energy on line 5 is ignored. 


The x and y source position 

parameters given in Table I-A are perpendicular to the axis of the detector shown in Figure 
1. z is along the axis. The program is now capable of modelling rectangular or circular 
sources, which results in a more accurate simulation of a source that is close to the detector 
(or widespread). This can be particularly important in studying the effects of collimation 
on a spectrum. To model a point source, set the maximum and minimum radius values to 
zero and the program will run as before. 

Table I-B shows the input and output files required to run BSIMUL. The isotope data files 
contain the energy and emission probability of each gamma-ray emitted with greater than 
0.1% probability. Up to twenty isotopes may be used, and the number may easily be 
increased. The isotopes whose energies and emission probabilities are available are listed in 
Table II. In order to produce a realistic spectrum, the simulation must include elements that 
form a part of the natural background. Terrestrial background elements that form a 
prominent part of any spectrum are 4(> K, 2l4 Bi, 2l2 Pb, and 228 Ac. The ratios of the activities 



Figure 1: A sample detector and shield configuration. 
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Table I-A. Input Parameters to BSIMUL. 


Input file line 

Comment 

N 

Y/N to plot the path of the gamma-rays through the 
detector system 

counts 

Event units: counts, nCi, pCi or mCi 

100 

Number of events 

5676654 9612341 

Random number seeds 

1000. 

Energy (keV) - used if no isotopes are specified 

Source Position Parameters 

39.68 

Source position on the z axis 

C 

Source symmetry: Rectangular or Circular (R or C) 

0. 2. 

Min, max source x dimension if rectangular, Radius if 

circular 


Min, max source y dimension (unused if circular) 

Energy Calibration Parameters 

1500. 

Maximum energy 

20. 

Minimum energy 

4096 

Maximum channel number 

50. 

Lower limit for detection in the shield 

5000. 

Lowest energy for Bremsstrahlung 

70. 

Effective charge for Bremsstrahlung 


of these elements relative to the activity of 2l4 Bi have been found to be: 

Isotope Activity relative to 214 Bi 

40 K 4.748 

2l4 Bi 1.0 

2l2 Pb 0.5227 

228 Ac 0.4998 


The simulation program uses the relative activities in calculating the probability of a 
Includes Constellation Technology Competition Sensitive Proprietary Data 
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background gamma-ray. 


The first output file will contain the spectrum of events occurring in the detector, the 
second will contain the spectrum of events occurring in the detector and not vetoed by the 
shield. The third and fourth files are used only when a collimator is specified. They 
contain spectra of events not vetoed by the collimator, and events vetoed neither by the 
collimator nor by the shield. 

Table I-C contains the detector configuration and materials of the Near Earth Asteroid 
Rendezvous (NEAR) detector. The detector part is specified with a D in the second 
column; the shield and collimator are specified with an S and a C, respectively. The next 
variable is the material. Except in the case of BGO, this is written using the chemical 
formula. While the collimator material and dimensions are entered here, the hole pattern 
must be entered in the program itself. 


Table I-B. Input Isotope and Output Spectrum Files. 


Input File Line 

Comment 

2 

Number of isotopes 

bi-214 0.5 

Isotope and percent probability of emission from that 

isotope 

ac-228 0.5 


0 . 

Percentage of background isotopes 

testbg2.sp 

File name of spectrum of all counts in detector 

testbgs2.sp 

File name of spectrum of all counts not vetoed by the 

shield 

temp5.sp 

File name of spectrum of all counts not vetoed by the 
collimator 

temp6.sp 

File name of spectrum of all counts not vetoed by the 
collimator or shield 
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Table I-C. Input Detector Material and Configuration Parameters. These Are the 
Parameters Used to Model the NEAR Detector. 


part 

material 

inner 

radius 

at 

bottom 

outer 

radius 

at 

bottom 

bottom 

inner 
radius 
at top 

outer 
radius 
at top 

top 

D 

Nal 

0.0 

1.26 

3.93 

0.0 

1.26 

12.36 


Al 

1.26 

1.57 

3.10 

1.26 

1.57 

12.36 


A1 

0.0 

1.26 

12.36 

0.0 

1.26 

12.67 


Al 

0.0 

1.26 

3.62 

0.0 

1.26 

3.93 

S 

BGO 

1.91 

4.45 

3.10 

1.91 

4.45 

10.68 

S 

BGO 

1.91 

4.45 

10.68 

WEm 

4.17 

14.68 

S 

BGO 

0.0 

4.45 

0.0 

0.0 

4.45 

3.10 
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Table II. Isotopes with Gamma-Rays Available to the Program. 



“Co 

85 Se 

85 Sr 

88 y 

95 Zr 

97 Zr 

89 Rb 

^Nb 

9, Nb 

95m Nb 

"Mo 

l03 Ru 

K)6r u 

l05 Tc 

108 a r pQ 

108b Tc 

“ 9 Cd 

ll9 In 

l27a In 

127b In 

l26 Sn 

l27 Sn 

l30 Sn 

l32 Sn 

132 Te 

l33 Te 

133m Te 

1 29 j 

,3, I 

132 j 

l33 Xe 

,34 Cs 

l36 Cs 

l37 Cs 

l40 Ba 

l42 Ba 

l43 Cs 

l46 La 

151 Sm 

2l2 Bi 

2l4 Bi 

ll9 Cd 

192 Ir 

208^1 

2l0 Pb 

21l Pb 

2l2 Pb 

2'4p b 

2l9 Rn 

226 Rn 

226 Ra 

227 Th 

228 Th 

228 Ac 

23 ’Pa 

233 Pa 

235 u 

238 U 

239 u 

239 Np 

24l Am 

243 Am 

244 Arn 

244m Am 

252 Cf 
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Absorption Coefficients 


The simulation program uses linear absorption coefficients in calculating interaction 
probabilities. A commercial program, PhotCoef, provides the linear photoabsorption, 
Compton, and pair production coefficients for each element in a compound. BSIMUL can 
use this to calculate absorption coefficients of any compound as follows: 


M; 

P 


N U, 

M P.. ' 


3 


where p s may represent the coefficient for photoabsorption, Compton absorption, or pair 
production of element i. N is the number of elements in the compound. Wj is the fraction 
by weight of the element i in the absorber, Pj is the density of the element, and p is the 
density of the compound. 

The absorption coefficients of compounds calculated by BSIMUL were found to be the 
same as the coefficients calculated with PhotCoef to five significant figures. The absorption 
constants for each material in the library, shown in Table III, cover an energy range of 10 
keV to 10 MeV in intervals of 10 keV up to 100 keV, in intervals of 100 keV up to 1 
MeV, and in intervals of 1 MeV up to 10 MeV. The library will be expanded to contain all 
of the elements, so that any material can be modeled. 


Table III. Materials with Linear Absorption Coefficients Calculated Using PhotCoef. 


Elements 

Li 

Be 

B 

C 

N 

O 

F 

Ne 

Na 

A1 

Si 

P 

S 

Cl 

K 

Cu 

Ge 

Cd 

Te 

Xe 

Hg 

Pb 

Bi 


Compounds 

Nal 

Csl 

... 

Hgl 2 

CdTe 

BGO 



Includes Constellation Technology Competition Sensitive Proprietary Data 


8 































Detector Efficiency 


The relative efficiency of a detector configuration can now be calculated from BSIMUL 
output. At the end of a simulation, the program now calculates the peak strengths by 
subtracting the average of twenty background points from the peak count value. The error 
is set equal to the square root of the average background value. Three files are produced. 

File Contents 

fn.xda Energy calibration constants 

fn.pk Peak position and error, FWHM and error, strength and error, 

and peak type 

fn.ca Peak position and error, energy and error, peak strength and 

error, emission probability and error 

These files are used as the input files to the calibration program XCALIBER, which 
calculates a fit to the efficiency using cubic splines. The efficiency curve produced in this 
way may be considered to be an upper limit to the actual detector efficiency, since the 
BSIMUL peaks are one channel wide, and no strength is lost in the long tail that is 
primarily due to nonideal detector response and some Compton scattering. 

One warning must be issued here. Both ‘ l4 Bi and 228 Ac emit a 904 keV gamma-ray. When 
214 Bi and 228 Ac were first used as the efficiency calibration elements, the whole strength of 
that peak was attributed to both elements, and the calibration results were wildly wrong. In 
writing the calibration file the program now omits peaks whose energies are close enough 
that their photopeaks are placed in the same channel. 


Calculations 


NEAR Nal Detector 

The configuration of the NEAR detector is shown in Figure 2. Several simulations have 
been performed with this configuration. The spectra were convoluted with a Nal peak 
shape. The spectra shown in Figure 3 are the BSIMUL output and convolution of a 214 Bi 
and 228 Ac spectrum. The output from BSIMUL was used to calculate the efficiency of the 
detector as a function of energy. The shield does lower the Compton background 
significantly, as is shown in Figure 5. The efficiency, shown in Figure 6, is the same for 
both the nonvetoed spectrum and the spectrum of all counts. A fit to the convoluted 
spectrum of nonvetoed counts, shown in Figure 7, produced the peaks listed in Table IV. 
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BGO Shield 


Figure 2. The NEAR Detector Configuration 


A series of runs were performed to study the count rate as a function of the source position. 
Figures 7-9 show three of the spectra, which were calculated with Co and Cs sources. The 
peak and background count values in the 661 keV peak as a function of source position are 
listed in Table V. The sharp reduction in the peak count rate for all the energies that 
occurs at 1 1.5 degrees is due to the source being moved into the shadow of the shield. The 
increase in the peak count rate as the source is moved away from the vertical axis from 40 
to 80 degrees has been seen experimentally in a germanium detector. As Figure 9 shows, 
the shield is quite effective in preventing gamma-rays from entering through the side of the 
detector. 
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Energy (keV) 

Figure 3: NEAR simulation. The spectrum of non-vetoed counts 
before and after convolution. Sources: Bi 214 and Ac 228 
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Energy (keV) 

Fiqure 5: NEAR spectrum. Fit to the spectrum of non-vetoed 
counts. Sources: Bi 214 and Ac 228. 



Efficiency 



Figure 6: Simulated NEAR Efficiency. 
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Table IV. Peaks Found in the Convoluted NEAR Spectrum that was Simulated with 

2,4 Bi and 228 Ac Isotopes. 


Analys 

is Results 


Isotope Data [Ref 2] 

Energy (keV) 

Peak strength 

Energy 

Emission 

Isotope 


(counts) 

(keV) 

Probability 






(%) 


59.82 +/- .5 

1398.12 

+/- 69.41 

57.81 

.525 

228 Ac 

92.99 +/- 1 . 

29922.5 

+/- 245.6 

89.955 

3.4 

228 Ac 




93.35 

5.6 

228 Ac 

127.83 +/- .5 

5929.66 

+/- 182.2 

129.03 

2.9 

228 Ac 

153.06 +/- 1.2 

1809.79 

+/- 141.4 

153.89 

.84 

228 Ac 

207.17 +/- .28 

7599.96 

+/- 124.9 

209.3 1 

4.1 

228 Ac 

268.04 +/- .4 

5016.51 

+/- 113.4 

270.26 

3.8 

228 Ac 

333.92 +/- .16 

18242.9 

+/- 168.8 

328.07 

3.5 

228 Ac 




338.42 

12.4 

228 Ac 

409.88 +/- 1 . 

1914.51 

+/- 92.24 

409.62 

2.2 

228 Ac 

461.29 +/- .7 

3960.20 

+/- 101.6 

463.1 

4.6 

228 Ac 

594.92 +/- 4. 

9288.07 

+/- 4570. 




609.68 +/- 2. 

21103.3 

+/- 4595. 

609.3 1 1 

46.1 

214 Bi 

751.88 +/- 5.6 

1429.28 

+/- 422.4 

755.28 

1.32 

228 Ac 




768.35 

4.88 

214 Bi 

783.20 +/- 2.4 

3972.13 

+/- 426.4 

794.79 

4.6 

228 Ac 

910.28 +/- .5 

12872.1 

+/- 221.6 

911.16 

29. 

228 Ac 

963.81 +/- .8 

8690.42 

+/- 211.5 

964.64 

5.8 

228 Ac 

1116.44 +/- .6 

6118.90 

+/- 108.0 

1120.27 

15. 

214 Bi 

1188.78 +/- 7. 

838.668 

+/- 83.50 

1 155.18 

1.69 

214 Bi 

1246.04 +/- 1.8 

2196.66 

+/- 108.4 

1238.11 

5.9 

2l4 Bi 

1373.92 +/- 1.2 

1634.47 

+/- 54.26 

1377.66 

4.02 

214 Bi 
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lotion. The spectrum of all counts and 
i— vetoed counts Sources: Co 60 and Cs 137 
detector axis. 
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Energy (keV) 

jlation. The spectrum of all counts and 
n— vetoed counts Sources: Co 60 and Cs 137 
8, theta = 1 6.6 
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Energy (keV) 

tion. The spectrum of all counts and 
vetoed counts Sources: Co 60 and Cs 137 
, theta = 60. 
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Figure 10: The peak count rates as a function of the angle between the source 
position and detector axis. 
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Table V. Source Position Dependance of the 661 keV Peak and Background Count 

Rate. 


X 

Theta 

(degrees) 

Counts in 661 
keV peak 

Counts in 
background 
of nonvetoed 
spectrum at 
661 keV 

Counts in 
background 
of spectrum 
of all counts 
at 661 keV. 

0. 

0. 

2870 

13 

91 

4. 

1 1.5 

1774 

9 

49 

8. 

16.6 

4422 

25 

144 

12. 

25 

1980 

16 

81 

16. 

33.4 

683 

5 

49 

20. 

42 

202 

3 

36 

24. 

50.7 

247 

5 

36 

28. 

60 

259 

3 

40 

32. 

70 

418 

8 

49 

36. 

82 

536 

5 

45 

40. 

107 

532 

9 

49 


Ge Detectors 

The work discussed in this section was performed pursuant to the requirements of contract 
F08650-94-C-0006, for the Advanced Research Projects Agency (ARPA) and the U.S. Air 
Force. It is included here as an illustration of the capabilities of the program. The work 
also proved to be useful as a full test of several of the modifications that have been made 
under the NASA contract. 

The germanium detector response was modeled with BS1MUL for various detector/shield 
geometries, and source/background conditions. Figure 10 shows an overlay of data taken 
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with a 30% Ge detector and a BGO shield. l37 Cs and '’"Co sources of roughly equal 
strength were used. A 1% background was assumed, and the shield shape was 
approximated, since the actual shield does not have cylindrical symmetry. The simulated 
spectrum was not given an accurate peak shape, since the interest was in the peak strengths 
only. 

The simulations were done using the dimensions for a 100% Ge crystal unless otherwise 
specified; and a 2 mm thick aluminum case around the crystal was assumed. A minimum 
of one million gamma-rays were used for each spectrum, a study of the detector 
performance as a function of shield thickness was made. Figure 1 1 shows the peak-to- 
Compton ratios and Figures-of-Merit obtained using Csl and BGO shields. The background 
was 10 times the source strength. The shield extended 6 cm in front of the detector, so that 
the acceptance angle was about 69". It was found that the optimum acceptance angle varies 
with the source strength and background count rate: if the shield extends far in front of the 
detector, the gamma-rays will Compton scatter from the shield to the detector. 

Simulations were run to determine the difference in performance between N and P type Ge 
detectors as a function of detector efficiency. The detector and shield mass were varied to 
keep the total mass constant, so as the detector efficiency was increased, the shield 
thickness was decreased. The Figure-of-Merit (FOM) 

e(£.) 

W = v 4 

[R(E )e 2 (E )i 


is useful in evaluating the suitability of a detector for a particular task. e(E) is the 
efficiency at energy E, e c2 (E) is the Compton efficiency due to a second gamma-ray, and R 
is the resolution. A plot of the FOM relative to the FOM obtained using and un-shielded 
100% N-type detector vs. the detector efficiency for both N and P type detectors is shown 
in Figure 12 for two low energy gamma-rays. At 129 keV the advantage of N over P type 
is evident, while at 375 keV there appears to be no significant difference. 

A collimator with a 10° angle of acceptance was modeled using a range of gamma-ray 
energies. The peak count rate as a function of source position is shown in Figure 16. 
Because a point source was used, there is a reduction in peak counts occurring below 5°. 
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Figure 11. Comparison of BSIMUL data and experimental 



Figure 12 Performance vs. shield thickness, 1000% background. 69 deg acceptance angle 
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Progress on Related Work 


ROBFIT 

The work on the spectral analysis package ROBFIT was required through contract F08650- 
94-C-0006, for the Advanced Research Projects Agency and the U.S. Air Force. A new 
program, UNKFIT, has been added to ROBFIT. Its purpose is to fit spectra taken with 
germanium detectors and produce a list of the isotopes present and their relative 
concentrations. UNKFIT is now in the integration and debugging phase. During the 
course of writing UNKFIT, we were able to realize an enormous simplification in the 
automation of ROBFIT. This simplification was accomplished by including the functions 
previously done with a sequence of modules used to tit the background, fit the peaks, and 
then identify the isotopes. In a further simplification, UNKFIT now uses detector-specific 
characteristics. Only one menu file is now required. 

These changes are expected to result in greater reliability and have been shown to reduce 
the total analysis time of a spectrum taken for 1 5 minutes to less than 5 minutes. In its 
present state, UNKFIT detects only the naturally-occuring background isotopes. Initially, 
these were identified at about the 5-sigma level in all but one of ten Monte-Carlo-generated 
test spectra. Current work includes debugging the UNKFIT module and adding isotopes of 
possible interest to the library. 

The test spectrum ztcase9.sp, shown in Figure 15, was a background spectrum taken for 72 
hours in a room containing a depleted fuel rod that was fully shielded by lead bricks. The 
spectrum has been fully analyzed and the isotopes and their relative activities are listed in 
Table VI. The program is now able to report the uncertainty with which an isotope is 
known to be missing from the spectrum as well as to be present. For instance, in the case 
of 239 U the results can be interpreted as evidence that this isotope is missing. 


Table VI. Isotopes Found in the Background Spectrum ztcase9.sp. 


Isotope 

Activity relative to 
2i4 Bi 

Error in relative 
activity 

2l4 Bi 

1. 

.00324 

228 Ac 

.4998 

.0112 

252 Cf 

.3119 

1.17 

244 Am 

.00258 

.0004 

243 A m 

.001224 

.00116 

24l Am 

.04849 

.0103 
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239 u 

.002484 

.0230 

239 Np 

.01564 

.00256 

235 U 

.0499 

.00537 

234 Th 

2.133 

.173 

233 Pa 

.002646 

.000964 

23 'Pa 

.02171 

.00717 

228 Th 

2.335 

.305 

227 Th 

.008366 

.00640 

226 Ra 

1.573 

.173 

2l9 Rn 

.05646 

.00654 j 

2> 4 p b 

.9114 

.0261 

212 Pb 

.5227 

.0277 

212 Bi 

.5849 

.0113 

21 'Pb 

.06563 

.00839 

2l0 Pb 

1.042 

.197 

ZOS'yj 

.1919 

.00484 

151 Sm 

269 

95.1 

146 La 

.009242 

.00102 

142 Ba 

.02338 

.00215 

137 Cs 

.008366 

.000386 

134 Pa 

1.329 

.0578 

l33 Te 

.007349 

.00165 

133m Te 

.002175 

.000306 

188 Xe 

.04282 

.00720 

130 Sn 

.0004989 

.000538 

l32 Sn 

.002874 

.000482 

132 I 

.01477 

.00293 
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] 29 j 

.1457 

.0328 

l27 In 

.003770 

.000780 

l27 Sn 

.005990 

.00112 

l26 Sn 

.06098 

.00894 

119 Cd 

.007802 

.000701 

ll9 In 

.07988 

.0183 

I08 Tc 

.01746 

.00160 

106 RU 

.5493 

.514 

105 Tc 

.04410 

.00817 

9, Zr 

.003118 

.000988 

95 Nb 

.03256 

.00391 

94 Nb 

.0009307 

.000210 

89 Rb 

.004132 

.00144 

88 y 

.0001841 

.000166 

85 Se 

.01185 

.00127 

85 Sr 

.001052 

.000394 

“Co 

.001630 

.000124 

40 k 

4.748 

.0361 
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taken with a 90% Ge detector over 




Plans 


Gamma-ray Spectrometer Modeling Studies 

During the coming year, the modifications to BSIMUL will be completed. A graphical 
interface that allows the operator to draw the detector would be a great improvement over 
the current menu system. Popup menus would allow the operator to choose the material of 
each part and the exact dimensions would be entered through a text window. New 
subroutines will be written that will allow the operator to model a detector and shield of 
any configuration. Since removing the symmetry requirements will slow the calculation, 
the subroutines that currently take advantage of cylindrical or rectangular symmetry will be 
used where possible. The convolution program will be integrated into BSIMUL, so that a 
spectrum will be given an accurate peak shape automatically. The result of these changes 
will be a program that is much easier to use than the current version. 

BSIMUL has been tested and works well at energies ranging from 50 keV to 1500 keV. 
Testing of the x-ray and high energy gamma-ray regions will be performed. The low 
energy range will be extended from 10 keV to 1 keV so that both x-ray fluorescence and 
gamma-ray systems can be simulated. The program can be used to model the behavior of 
any desired detector system. 

The x-ray simulations will be verified with spectra obtained using a Hgl 2 detector. The 
detector will be used in conjunction with active and passive shielding materials to provide a 
thorough test of the program at x-ray energies. 
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